{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "np.set_printoptions(suppress=True)\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we will construct binary and multiclass logistic regression models. We will try our binary model on the {doc}`breast cancer </content/appendix/data>` dataset and the multiclass model on the {doc}`wine </content/appendix/data>` dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import data\n",
    "cancer = datasets.load_breast_cancer()\n",
    "X = cancer['data']\n",
    "y = cancer['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first define some helper functions: the logistic function and a standardization function, equivalent to `scikit-learn`'s `StandardScaler`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic(z):\n",
    "    return (1 + np.exp(-z))**(-1)\n",
    "\n",
    "def standard_scaler(X):\n",
    "    mean = X.mean(0)\n",
    "    sd = X.std(0)\n",
    "    return (X - mean)/sd "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The binary logistic regression class is defined below. First, it (optionally) standardizes and adds an intercept term. Then it estimates $\\boldsymbol{\\beta}$ with gradient descent, using the gradient of the negative log-likelihood derived in the {doc}`concept section </content/c3/s1/logistic_regression>`,\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}(\\boldsymbol{\\beta})}{\\partial \\boldsymbol{\\beta}} = \\frac{\\partial - \\log L (\\boldsymbol{\\beta})}{\\partial \\boldsymbol{\\beta}} = -\\mathbf{X}^\\top (\\mathbf{y} - \\mathbf{p}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BinaryLogisticRegression:\n",
    "    \n",
    "    def fit(self, X, y, n_iter, lr, standardize = True, has_intercept = False):\n",
    "        \n",
    "        ### Record Info ###\n",
    "        if standardize:\n",
    "            X = standard_scaler(X) \n",
    "        if not has_intercept:\n",
    "            ones = np.ones(X.shape[0]).reshape(-1, 1)\n",
    "            X = np.concatenate((ones, X), axis = 1)\n",
    "        self.X = X\n",
    "        self.N, self.D = X.shape\n",
    "        self.y = y\n",
    "        self.n_iter = n_iter\n",
    "        self.lr = lr\n",
    "\n",
    "        ### Calculate Beta ###\n",
    "        beta = np.random.randn(self.D) \n",
    "        for i in range(n_iter):\n",
    "            p = logistic(np.dot(self.X, beta)) # vector of probabilities \n",
    "            gradient = -np.dot(self.X.T, (self.y-p)) # gradient\n",
    "            beta -= self.lr*gradient \n",
    "            \n",
    "        ### Return Values ###\n",
    "        self.beta = beta\n",
    "        self.p = logistic(np.dot(self.X, self.beta)) \n",
    "        self.yhat = self.p.round()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following instantiates and fits our logistic regression model, then assesses the in-sample accuracy. Note here that we predict observations to be from class 1 if we estimate $P(Y_n = 1)$ to be above 0.5, though this is not required. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In-sample accuracy: 0.9894551845342706\n"
     ]
    }
   ],
   "source": [
    "binary_model = BinaryLogisticRegression()\n",
    "binary_model.fit(X, y, n_iter = 10**4, lr = 0.0001)\n",
    "print('In-sample accuracy: '  + str(np.mean(binary_model.yhat == binary_model.y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the graph below shows a distribution of the *estimated* $P(Y_n = 1)$ based on each observation's *true* class. This demonstrates that our model is quite confident of its predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEiCAYAAAACg5K6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debhcVZnv8e+PJBAZGhIyGBI0jF5ChEOIgG2DDAoRbzM0qNAtJggEBWxo29vaeL2itIINgnob0SB0Ai0y2UDkMjREIDiABEgggEKEACeJIQQSoUNChvf+sdaBolKnqs6pqjNsfp/nqaeq1l5773dX7Xpr1dpr71JEYGZmxbJJbwdgZmbN5+RuZlZATu5mZgXk5G5mVkBO7mZmBeTkbmZWQE7uZmYF5ORuZlZATu5lJE2RFJ3cVnRhOUdJ+mKF8nMk9fiZY53F04Tl3iPpnjrqnVL2Wq6WNF/S5Ap1/6+kX5Q8v1fSYkmbVqjb8X6d3PDGdIGkMTnO30palWMYW6HeP0h6VFJdn7WO/UPSwBbE3Nl+XXpb2Oz1NkLSByVdl9//NyQtl3SnpMmSBuQ6vfKZ6uuc3Dv3CeCDZbePdGH+o4BKyfQneVk9rbN4ekobsJq3XsujgT8D0yUd1FFJ0k7AqcA3SuY9HxgFfKp0gZL2B34MfDciftLS6De2M/BJ4BXgvir1fgSMADb6EusF5fvzn4A7ysqO7rXoykg6C/g1MBT4Munz91ngKeBS4H/2XnR9X9NbBwUyNyIWNHuhEdEOtDd7uf1AG/BERNzfUZBbiU8AhwN35+KzgHkRMaejXkTcJulR4B+Aq/K8OwL/CfwX8E89sQFlZkfEyBzLycChlSpFxOuSrgS+BPx7D8ZXKZb7S59LWgO8VF7eGUmbRcSalgS38boOAC4C/i0i/r5s8s2SLgK26IlY+iu33LtJ0q6SbpT0Yu5ieF7S9ZIGSppOaqmNLv+5W+knZMlP8f8h6Q5J/52Xd2KefoKk30t6TdLduXXbMe/Okq6S9Kyk1yU9I+lSSUNK6nQaT56+p6SZkl7Jy/h1bhWXb/NxOY41kh6XVFcrT5KAPYDHyib9Od+/K9fbDPg0cHWFxfwrsJekAyVtDdwCLAGOj4gN9cTRTF1c5zXAOEl/2YV5dsvv9SpJSyR9s6NrR9Kx+T3cs3ym3E322y6sp6KSfXJ83idfA67L06ZX6r5RhS66evetCr4CvEwnX9wR8ceIeLRK/PV8Ljr9DHelTl/V5wPsRQMqvIEbSj7UtwArgM8DLwGjSS3QTYBzgeHAB4Ajcv16WjzXA5cBFwKnAVdI2gU4kLSzDwK+T0p+++Z5tiP9EjiL1EWwI3A2cCtvdf90Go+kCaRuhUeAU4BVwOeAuyT9ZUQ8lOt9JK/3/wH/mJf3/RzTH2ps1y7AlkD5h/HD+f6hfL8fsA2VuzmuBb4F/C/SfjsU2CciXqux7oryF86AOqpGRKzvzjpKzCV9kU0CflPnPDcBVwDnAYcBXwM2AOfkaYtJ3Vendcwg6X2k1/TEBuMtdTNwOfCdvP661btvVZhvAGmfvykiVncz7no+F9U+w3ShTt8UEb6V3IApQHRyuyXXGZafH1FlOdOB9grl56SXfeMy4DMlZUOAdcBy4C9Kyv8+131vJ+sdCPxVrrNXHfHMAp4ENi0pG5DLbiop+zWpC2WTkrJ983ruqfGafjLXm5TjGwL8DbA0r2dwrvdlUgLZtJPlfCEv53Vgvwbf5wOrvM+lt6rblpd1cq47tkqd+4D/qmNZHfvCV8rKLwNeBbYpqbcS2KKkzkWkRPauOl+DhcB/1IjjzE727YUVyu8pfb3q3bcqLGdkXvd5dW7HOZR9pmp9LqjvM1yzTl++ueXeuaPZuG+8Y7TMcuAZ4HxJI0k79NNNWOdtHQ8i4hVJLwKPRMSfS+r8Pt9vDzynNILkS8BngPcCg0vqvo/UaqpI0rtILb1vAxvKfqncBfxdrjeA1Oo/P0q6IyLigUo/zyvYq3z7gLXAjaTk0dE62w74c0S80clybgF+AFwYdfYTV/EQaZtqebXB9XRYBuzahfrXlT2/hvQlMh74FTAN+CpwPPATSYNJXW9XRsTrjYf7phu7M1O9+1ar1PG5mEvtz3CrPuc9wsm9c/OjkwOqERGSPkpqMZwHbCvpWeCCiLi0gXW+Uvb8jU7K4K2d9TxSi/abpJ/8rwJjSAcbB1PdUFJL6mv5tpHczzuM1P2ytEKVSmXl2khflEfzVsv72QpJaDDVu692y/e/q2OdtbxG+oDX0qwhdq+Tjy3Uqfx17Xg+GiAiFku6mdTN8RPS6K6hpNFDzbSkm/PVtW9F5WMXy0mv13u7uW6o8bmo5zPcws95j3By76aIeAb4TO673RM4A/ihpIURcVv1uZvqOFJr7V86CiRtWee8K0jdIJcAV1aqEBEbJL1EammPrFBlJPBcjfW0AfdHyQiYTiwnddlUWw508mtE0rdJX0TDSL8WFgGTonK//Id5a4RONfeSunAaNZTUZ1uvkaRWY+lzSNvU4YfALEl7k/rf74uIJxqKcmOVvtxWAxudcwBsS3oPoc59q5PydfnA7EfV/RE6NT8X9XyG+9DnvMv6/kGBPi6Subw1hnx8vl9D11pq3bU5KfGWqnRAbaN4IuK/SX3BewIPR8Sc8luutx54EDhWJSfjSNoXGFstuPxz9t1U6R4q8XtgkKQxnUzfE1geaThpJRNI4+E/TTqAFsABndTt6JapdTu1jrjrsQO1DzyX+mTZ8+NIvzbmdxRExC9J/dcXAR8ijanvCc8BIyUN6yhQGsH1vpLY6tq3qjif9GVxQaWJknaQtEeV+ev9XFT7DHepTl/jlnvn2kp33hJzgHGkkSLXAgtIPz+nkA6A/jLXewIYKunzeZ7VEVE+FLAZbgcmS3osx/I3QKUhd53F80VgNnCHpMtJP8OHkRLlgIj4Sp7/66Qx5TdJ+jFptMw3SCfCVNPR315Pcp+d7/eh8rkAbVTvSpkAfDgiVsGb/a4vV6oYEa+SXoduk3Rsfrh3vv+YpGXAsoi4t6TeNqT+9gu7sPhT8hfpg6TRMicD50RE+VnSPyLtiy8BP+/6VnTL9aQRWD9VGm8+DPhnNv5lUu++tZGImK10RvVFknYjHcR9nvTL7hDS6/G3bDwCq0PVz0X+Yqj6Ga6nTp/W20d0+9qN6qNlgrRzjgBmkM6UW0VKIPcCh5UsZwvgZ6Q+8yCPLqD6aJmBZeULKRvNwFujPD4Sbx3Rvyav5xXgp6QWZwBTasWTp+2Wl/EiqYXfDswEDi9b9/Gk1uca4HFSH/o9VBlRQhrCGcCYOl//B4B/r1C+ObCedDC10nzvISXVjueDSC3dukaOdHNfqWuEDeng4Wpg2zqW2bEvjCd1G71O+gI9l5KRSiX1R+X6F3Qj/o32r1r7ZMn0o0i/Il4H5pFO4tpoX6h336oS41+SvkyWkFriL5MaGZ/ueD2o/Jmq+rmgvs9wzTp9+aa8EWZ9gqQppNbSqMgt8DrnOwqYGhGH5+cTgOkRUe2ne4+QdBvpTNATWrDsU0gHUXeNFpxRbf2X+9ytr7mKdNDwtFoVy+zN27tZJtJgt0szSGoDDuLt18ppxnLHSfrrvNybnNitnFvu1udI2g+YEBE/7O1YGiVpEjAkIn7W5OXeQ+qy+A3wtxGxuJnLt/7Pyd3MrIDcLWNmVkBO7mZmBdQnxrlPmjQpbr/99t4Ow8ysv1FnE/pEy/2ll7pyVraZmdXSJ5K7mZk1l5O7mVkBObmbmRWQk7uZWQE5uZuZFZCTu5lZATm5m5kVkJO7mVkB9YkzVM3M+pNp0+Y1bVlTp+7ZtGWVcsvdzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCqhmcpc0WNLvJM2T9Likb+TyHSQ9IOlpSddK2jSXb5afL8jTx7Z2E8zMrFw9Lfc1wMERsSfQBkyStB/wHeDiiNgFeAU4Kdc/CXglInYGLs71zMysB9VM7pG8lp8OyrcADgZuyOUzgKPy4yPzc/L0QySpaRGbmVlNdfW5SxogaS7wInAn8EdgRUSsy1XagdH58WjgBYA8fSWwbTODNjOz6upK7hGxPiLagDHAPsBularl+0qt9CgvkDRV0hxJc5YtW1ZvvGZmVocujZaJiBXAPcB+wDaSOq4HPwZYnB+3A9sD5OlbAy9XWNa0iJgYEROHDx/evejNzKyimn/WIWk4sDYiVkh6F/AR0kHSu4FjgWuAycDNeZaZ+flv8/RfRsRGLXcz69vWrl1Le3s7q1ev7u1Q+pxdd4VnngnWreu7hxPr+SemUcAMSQNILf3rIuIWSU8A10j6F+AR4PJc/3LgKkkLSC3241oQt5m1WHt7O1tttRVjx47FYyLeEhFsuuli4CWeeqq3o+lczeQeEY8Ce1Uof4bU/15evhr4RFOiM7Nes3r1aif2CiSx1VbbsPnmL/V2KFX5DFUz65QTe2X94XVxcjezPutPf/oTxx13HDvttBPjxo3j8MMP56mnnmLhwoWMHz++Jetcs2YNn/rUp9h5553Zd999WbhwYUvW02r19LmbmTH78eYeWD1g98FVp0cERx99NJMnT+aaa64BYO7cuSxdupTtt9++qbGUuvzyyxkyZAgLFizgmmuu4ctf/jLXXntty9bXKm65m1mfdPfddzNo0CA+97nPvVnW1tbG/vvv/7Z6CxcuZP/992fChAlMmDCB3/zmNwAsWbKEAw44gLa2NsaPH899993H+vXrmTJlCuPHj+f9738/F1988Ubrvfnmm5k8eTIAxx57LLNmzaI/Dvhzy93M+qT58+ez995716w3YsQI7rzzTgYPHszTTz/N8ccfz5w5c7j66qs57LDD+OpXv8r69etZtWoVc+fOZdGiRcyfPx+AFStWbLS8RYsWvfnLYODAgWy99dYsX76cYcOGNXcDW8zJ3cz6tbVr13LGGWcwd+5cBgwYwFN5fOIHPvABPvvZz7J27VqOOuoo2tra2HHHHXnmmWf4whe+wMc//nEOPfTQjZZXqZXeHw6glnO3jJn1SbvvvjsPPfRQzXoXX3wxI0eOZN68ecyZM4c33ngDgAMOOIDZs2czevRoTjjhBK688kqGDBnCvHnzOPDAA7nkkks4+eSTN1remDFjeOGFFwBYt24dK1euZOjQoc3duB7g5G5mfdLBBx/MmjVruOyyy94se/DBB7n33nvfVm/lypWMGjWKTTbZhKuuuor169cD8NxzzzFixAhOOeUUTjrpJB5++GFeeuklNmzYwDHHHMO5557Lww8/vNF6jzjiCGbMSBe2veGGGzj44IP7Zcvd3TJm1idJ4sYbb+Sss87i/PPPZ/DgwYwdO5bvfe97b6t32mmnccwxx3D99ddz0EEHscUWWwBwzz33cMEFFzBo0CC23HJLrrzyShYtWsSJJ57Ihg0bADjvvPM2Wu9JJ53ECSecwM4778zQoUPfHKnT36gvHAWeOHFizJkzp7fDMLMSTz75JLvtVukCsLZs2SqeffZp5s5tfFlTp+7ZyOyd/qRwt4yZWQE5uZuZFZCTu5lZATm5m5kVkJO7mVkBObmbmRWQk7uZ9Vm9ccnf2bNnM2HCBAYOHMgNN9zQknX0BJ/EZGZ1mTZtXlOXV2t8d29d8vc973kP06dP58ILL2zZOnqCW+5m1if11iV/x44dyx577MEmm/Tv9OiWu5n1Sb11yd+icHI3s36t2Zf8LYr+/bvDzAqrty75WxRO7mbWJ/XWJX+LomZyl7S9pLslPSnpcUln5vJzJC2SNDffDi+Z558lLZD0B0mHtXIDzKyYOi75e+edd7LTTjux++67c84557Dddtu9rd5pp53GjBkz2G+//XjqqafedsnftrY29tprL37+859z5plnsmjRIg488EDa2tqYMmVKxUv+Pvjgg4wZM4brr7+eU089ld13371HtrfZal7yV9IoYFREPCxpK+Ah4Cjgk8BrEXFhWf1xwM+AfYDtgLuAXSNifWfr8CV/zfoeX/K3c4W45G9ELImIh/PjV4EngdFVZjkSuCYi1kTEs8ACUqI3M7Me0qU+d0ljgb2AB3LRGZIelXSFpCG5bDTwQsls7VT/MjAzsyarO7lL2hL4OXBWRPwZuBTYCWgDlgDf7ahaYfaN+n4kTZU0R9KcZcuWdTlwMzPrXF3JXdIgUmL/aUT8J0BELI2I9RGxAbiMt7pe2oHSc4PHAIvLlxkR0yJiYkRMHD58eCPbYGYt0hf+hrMv6g+vSz2jZQRcDjwZEReVlI8qqXY0MD8/ngkcJ2kzSTsAuwC/a17IZtYTBg8ezPLly/tFIutJEcGrr65g1arejqS6es5Q/RBwAvCYpI5jw2cDx0tqI3W5LAROBYiIxyVdBzwBrANOrzZSxsz6pjFjxtDe3o67TTe2dOkbPPNMUGWwSq+rmdwj4ldU3oJbq8zzLeBbDcRlZr1s0KBB7LDDDr0dRp90333z6MuJHXyGqplZITm5m5kVkJO7mVkBObmbmRWQk7uZWQE5uZuZFZCTu5lZATm5m5kVkJO7mVkBObmbmRWQk7uZWQE5uZuZFZCTu5lZATm5m5kVkJO7mVkBObmbmRWQk7uZWQE5uZuZFZCTu5lZATm5m5kVkJO7mVkBObmbmRWQk7uZWQHVTO6Stpd0t6QnJT0u6cxcPlTSnZKezvdDcrkk/UDSAkmPSprQ6o0wM7O3q6flvg74x4jYDdgPOF3SOOArwKyI2AWYlZ8DfAzYJd+mApc2PWozM6uqZnKPiCUR8XB+/CrwJDAaOBKYkavNAI7Kj48ErozkfmAbSaOaHrmZmXWqS33uksYCewEPACMjYgmkLwBgRK42GnihZLb2XGZmZj2k7uQuaUvg58BZEfHnalUrlEWF5U2VNEfSnGXLltUbhpmZ1aGu5C5pECmx/zQi/jMXL+3obsn3L+bydmD7ktnHAIvLlxkR0yJiYkRMHD58eHfjNzOzCuoZLSPgcuDJiLioZNJMYHJ+PBm4uaT8M3nUzH7Ayo7uGzMz6xkD66jzIeAE4DFJc3PZ2cD5wHWSTgKeBz6Rp90KHA4sAFYBJzY1YjMzq6lmco+IX1G5Hx3gkAr1Azi9wbjMzKwBPkPVzKyA6umW6dNmP766ofkP2H1wkyIxM+s73HI3MysgJ3czswJycjczKyAndzOzAnJyNzMrICd3M7MCcnI3MysgJ3czswJycjczKyAndzOzAnJyNzMrICd3M7MCcnI3MysgJ3czswJycjczKyAndzOzAnJyNzMrICd3M7MCcnI3MysgJ3czswJycjczK6CayV3SFZJelDS/pOwcSYskzc23w0um/bOkBZL+IOmwVgVuZmadq6flPh2YVKH84ohoy7dbASSNA44Dds/z/FDSgGYFa2Zm9amZ3CNiNvByncs7ErgmItZExLPAAmCfBuIzM7NuaKTP/QxJj+ZumyG5bDTwQkmd9lxmZmY9qLvJ/VJgJ6ANWAJ8N5erQt2otABJUyXNkTRn2bJl3QzDzMwq6VZyj4ilEbE+IjYAl/FW10s7sH1J1THA4k6WMS0iJkbExOHDh3cnDDMz60S3krukUSVPjwY6RtLMBI6TtJmkHYBdgN81FqKZmXXVwFoVJP0MOBAYJqkd+DpwoKQ2UpfLQuBUgIh4XNJ1wBPAOuD0iFjfmtDNzKwzNZN7RBxfofjyKvW/BXyrkaDMzKwxPkPVzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKqGZyl3SFpBclzS8pGyrpTklP5/shuVySfiBpgaRHJU1oZfBmZlZZPS336cCksrKvALMiYhdgVn4O8DFgl3ybClzanDDNzKwraib3iJgNvFxWfCQwIz+eARxVUn5lJPcD20ga1axgzcysPt3tcx8ZEUsA8v2IXD4aeKGkXnsu24ikqZLmSJqzbNmyboZhZmaVNPuAqiqURaWKETEtIiZGxMThw4c3OQwzs3e27ib3pR3dLfn+xVzeDmxfUm8MsLj74ZmZWXd0N7nPBCbnx5OBm0vKP5NHzewHrOzovjEzs54zsFYFST8DDgSGSWoHvg6cD1wn6STgeeATufqtwOHAAmAVcGILYjYzsxpqJveIOL6TSYdUqBvA6Y0GZWZmjfEZqmZmBeTkbmZWQE7uZmYF5ORuZlZATu5mZgXk5G5mVkBO7mZmBeTkbmZWQE7uZmYF5ORuZlZATu5mZgXk5G5mVkBO7mZmBeTkbmZWQE7uZmYF5ORuZlZATu5mZgXk5G5mVkBO7mZmBeTkbmZWQE7uZmYF5ORuZlZATu5mZgU0sJGZJS0EXgXWA+siYqKkocC1wFhgIfDJiHilsTDNzKwrmtFyPygi2iJiYn7+FWBWROwCzMrPzcysB7WiW+ZIYEZ+PAM4qgXrMDOzKhpN7gH8l6SHJE3NZSMjYglAvh9RaUZJUyXNkTRn2bJlDYZhZmalGupzBz4UEYsljQDulPT7emeMiGnANICJEydGg3GYmVmJhlruEbE4378I3AjsAyyVNAog37/YaJBmZtY13U7ukraQtFXHY+BQYD4wE5icq00Gbm40SDMz65pGumVGAjdK6ljO1RFxu6QHgesknQQ8D3yi8TDNzKwrup3cI+IZYM8K5cuBQxoJyszMGuMzVM3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgBr9mz0zs35h2rR5vR1Cj3LL3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKqGUnMUmaBHwfGAD8JCLOb9W6zKxxzTzJZ+rUPZu2rHfayUfN0pLkLmkAcAnwUaAdeFDSzIh4ohXr6y2zH1/da+s+YPfB3Z63kbgbWW9/11eTX1/khNz7WtVy3wdYEBHPAEi6BjgSKFRy7029+cXSCH+xNJeTqHVGEdH8hUrHApMi4uT8/ARg34g4o6TOVGBqfvo+4A/dXN0w4KUGwu2PvM3vDN7md4ZGtvmliJhUaUKrWu6qUPa2b5GImAZMa3hF0pyImNjocvoTb/M7g7f5naFV29yq0TLtwPYlz8cAi1u0LjMzK9Oq5P4gsIukHSRtChwHzGzRuszMrExLumUiYp2kM4A7SEMhr4iIx1uxLprQtdMPeZvfGbzN7wwt2eaWHFA1M7Pe5TNUzcwKyMndzKyA+k1ylzRJ0h8kLZD0lQrTN5N0bZ7+gKSxPR9lc9WxzV+U9ISkRyXNkvTe3oizmWptc0m9YyWFpH4/bK6ebZb0yfxePy7p6p6Osdnq2LffI+luSY/k/fvw3oizWSRdIelFSfM7mS5JP8ivx6OSJjS80ojo8zfSQdk/AjsCmwLzgHFldU4DfpQfHwdc29tx98A2HwRsnh9//p2wzbneVsBs4H5gYm/H3QPv8y7AI8CQ/HxEb8fdA9s8Dfh8fjwOWNjbcTe4zQcAE4D5nUw/HLiNdI7QfsADja6zv7Tc37ycQUS8AXRczqDUkcCM/PgG4BBJlU6m6i9qbnNE3B0Rq/LT+0nnE/Rn9bzPAOcC/wr0z2swvF0923wKcElEvAIQES/2cIzNVs82B/AX+fHW9PPzZCJiNvBylSpHAldGcj+wjaRRjayzvyT30cALJc/bc1nFOhGxDlgJbNsj0bVGPdtc6iTSN39/VnObJe0FbB8Rt/RkYC1Uz/u8K7CrpF9Luj9fcbU/q2ebzwE+LakduBX4Qs+E1mu6+nmvqWWX/G2ympczqLNOf1L39kj6NDAR+HBLI2q9qtssaRPgYmBKTwXUA+p5nweSumYOJP06u0/S+IhY0eLYWqWebT4emB4R35X0QeCqvM0bWh9er2h6/uovLfd6LmfwZh1JA0k/5ar9DOrr6rqEg6SPAF8FjoiINT0UW6vU2uatgPHAPZIWkvomZ/bzg6r17ts3R8TaiHiWdJG9XXoovlaoZ5tPAq4DiIjfAoNJF9gqqqZfsqW/JPd6LmcwE5icHx8L/DLykYp+quY25y6KH5MSe3/vh4Ua2xwRKyNiWESMjYixpOMMR0TEnN4Jtynq2bdvIh08R9IwUjfNMz0aZXPVs83PA4cASNqNlNyX9WiUPWsm8Jk8amY/YGVELGloib19FLkLR5sPB54iHWX/ai77JunDDenNvx5YAPwO2LG3Y+6Bbb4LWArMzbeZvR1zq7e5rO499PPRMnW+zwIuIv0fwmPAcb0dcw9s8zjg16SRNHOBQ3s75ga392fAEmAtqZV+EvA54HMl7/El+fV4rBn7tS8/YGZWQP2lW8bMzLrAyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNytx0iaLqlfXRNG0i2SpjdhOUMkLZW0UxPC6jGSbpD0xd6Ow7rOyd0qyok4Ktzur3P+eyT9W1nxmcCnmx9tXetu9Tp/WfIarZP0R0mnllQ5G7g1Iv6Y6w+Q9JSkSyss6zuSFrf6+vySDpA0U9KiHPeUCtW+AfxvSVu3MhZrPid3q+YuYFTZrdt/mhDp8gH99WJXtexFupLhKNJ1X24DLpW0l6TNgZOByzsqR8R64NvAFEkjO8olnUa6Nv/HI+K5Fse8JTCf9KX7eqUKEfEY6VIHLf9StuZycrdq1kTEn8puL8Obrb77Jb0maaXSv1+Nz9Omk65QeXpJa3ZsebdMbmFfKum7kl6WtEzSmUr/qnWJpBWSnpd0QmlQ+V987pP0Sp7vjnz9kU7XnadJ0j/lVvXrkh7LV9TsWO7mOcbXchfK2fW8SLmrZRvg1/k1epZ0zXmRkv7hwAbS6fSl/oN0SvpZeTlHAt8FjomIR+pZdyMi4taIODsibsjxdWYm6SqN1o84uVuX5atu3gz8CtgT2Bf4PrA+VzkT+C3w77zV4n9h4yUB8HfAq3kZ5/YYgRsAAARUSURBVAPfI10o6ynSZYxnAD+RtF3JPFvkevuQLoO7EvhFvghVtXX/C+maHqeTrl1yHvBjSR/P0y8EPgocQ7po1V6kf9CpZe98P6+krOOPU14E9gceirJrfUT634Hzgc9LOhS4Gjg1Iu6sY50ASDo7fxlVu+1f7/I68TtgH0nvanA51pN6+4I6vvXNGzAdWAe8Vnb7DjCUdK3pD1eZ/x7g3yos85ayOr8teS7Slf9mlpQNAt4Ajq2yri1IXyx/VWXdW5C6HvYvK/8e6c8gtgTWAH9XMm1LYAXpuuLVXqvvAItKnu9E+uJbQLqg3U3AjE7m3ZR0BcT1wNndeJ+GAjvXuL2rjuW8BkzpZNoe+f3eqbf3S9/qv/WXP+uw3jEbmFpWtiIiXs7dH3dImgXMAq6PiM5a59U82vEgIkLSi6Sr4nWUrZX0CjCioyx3g5xLau0PJ/0C3QR4T5X1jCMl2tsllbagBwELSQl5U1Krv2Pdr0l6jNr2Bt4t6TXS/4MKuBE4PiJW5xbv0kozRsQbku4C3h8R365jXeXzv0zr/7egoz/eLfd+xMndqlkVEQsqTYiIEyV9D5gEHAF8S9JREXFHF9extnzRnZSVdiH+AlgEnJrv15Euh7tplfV0zP/XpJZyeQxD6g95I3uRLsn7Y1IiXBJv/8egl2osfw/StenfRtIvSH+M/VHSHzkcEREPl9U5mzQSp5qPRcR9tTaiiqH5vsjXUy8cJ3frtoiYR+pn/o6k20h/ltKR3N8gtWKbStK2wG7A6RFxdy6bwNv35UrrfoLU7fLeiPhlheUuJyX5/ch/hCFpC9I/P/2xSjw7kJLfXZ19EZIS9JRO5h+U13FJhcnjSd1YH8yjaI4GHi6r8yPyPxZVsajG9FrGA4sjouKvD+ubnNytms0kvbusbD2pL/pU0iiKRcCOpNZn6ZjthaSDcGNJ/bnN6jp4hdQSPkXSC6Q/Eb6A1HrvdN0R8aqkC4ELJYnU5bQlKZlviIhpki4nfVEtI/3F2f+h9hdUx8HUav8GdUde7rYRsbxs2nhgM8qStqStgE0i4se5aABQPm9D3TKStiT1yUPu1pLURnq9Sn/d7A/c3p11WO/xaBmr5iOkoXqlt0eAVaS/erueNKplBvBT0oHFDheSWtBPkH7OV+sPr1vu7vgU6ctkPqnF+zVSq7zWur9GGov+JeBx4E7SyJhn8/QvAXeT+svvzsufXSOkvYFnKyTt0pgfI404Oa7C5Ak59sfLyseT/o6uw/tzPM00kfR+PkLqT/9GfvzNjgqSBpN+MVzW5HVbi/mfmMx6gKRJpOGi4yKdwFSr/inAuyPi3Pz8fuDInu4akXR6Xu+hPblea5xb7mY9ICJuJ/3KGFOrbvZ+8kii3I00opf6vNcCX+iF9VqD3HI3Mysgt9zNzArIyd3MrICc3M3MCsjJ3cysgJzczcwKyMndzKyAnNzNzArIyd3MrICc3M3MCuj/Axdr8yoyjgYMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.distplot(binary_model.p[binary_model.yhat == 0], kde = False, bins = 8, label = 'Class 0', color = 'cornflowerblue')\n",
    "sns.distplot(binary_model.p[binary_model.yhat == 1], kde = False, bins = 8, label = 'Class 1', color = 'darkblue')\n",
    "ax.legend(loc = 9, bbox_to_anchor = (0,0,1.59,.9))\n",
    "ax.set_xlabel(r'Estimated $P(Y_n = 1)$', size = 14)\n",
    "ax.set_title(r'Estimated $P(Y_n = 1)$ by True Class', size = 16)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiclass Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import data\n",
    "wine = datasets.load_wine()\n",
    "X = wine['data']\n",
    "y = wine['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before fitting our multiclass logistic regression model, let's again define some helper functions. The first (which we don't actually use) shows a simple implementation of the softmax function. The second applies the softmax function to each row of a matrix. An example of this is shown for the matrix \n",
    "\n",
    "$$\n",
    "\\mathbf{Z} = \\begin{bmatrix} 1 & 1 \\\\ 0 & 1 \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "The third function returns the $I$ matrix discussed in the concept section, whose $(n, k)^\\text{th}$ element is a 1 if the $n^\\text{th}$ observation belongs to the $k^\\text{th}$ class and a 0 otherwise. An example is shown for \n",
    "\n",
    "$$\n",
    "\\mathbf{y} = \\begin{bmatrix} 0 & 0 & 1 & 1 & 2 \\end{bmatrix}^\\top.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Softmax for Z:\n",
      " [[0.5  0.5 ]\n",
      " [0.27 0.73]]\n",
      "I matrix of [0,0,1,1,2]:\n",
      " [[1 0 0]\n",
      " [1 0 0]\n",
      " [0 1 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def softmax(z):\n",
    "    return np.exp(z)/(np.exp(z).sum())\n",
    "\n",
    "def softmax_byrow(Z):\n",
    "    return (np.exp(Z)/(np.exp(Z).sum(1)[:,None]))\n",
    "\n",
    "def make_I_matrix(y):\n",
    "    I = np.zeros(shape = (len(y), len(np.unique(y))), dtype = int)\n",
    "    for j, target in enumerate(np.unique(y)):\n",
    "        I[:,j] = (y == target)\n",
    "    return I\n",
    "\n",
    "\n",
    "Z_test = np.array([[1, 1],\n",
    "              [0,1]])\n",
    "print('Softmax for Z:\\n', softmax_byrow(Z_test).round(2))\n",
    "\n",
    "y_test = np.array([0,0,1,1,2])\n",
    "print('I matrix of [0,0,1,1,2]:\\n', make_I_matrix(y_test), end = '\\n\\n')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multiclass logistic regression model is constructed below. After standardizing and adding an intercept, we estimate $\\hat{\\mathbf{B}}$ through gradient descent. Again, we use the gradient discussed in the {doc}`concept section </content/c3/s1/logistic_regression>`,\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}(\\mathbf{B})}{\\partial \\mathbf{B}} = \\frac{\\partial -\\log L (\\mathbf{B})}{\\partial \\mathbf{B}} = \\mathbf{X}^\\top (\\mathbf{I} - \\mathbf{P}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MulticlassLogisticRegression:\n",
    "    \n",
    "    def fit(self, X, y, n_iter, lr, standardize = True, has_intercept = False):\n",
    "        \n",
    "        ### Record Info ###\n",
    "        if standardize:\n",
    "            X = standard_scaler(X) \n",
    "        if not has_intercept:\n",
    "            ones = np.ones(X.shape[0]).reshape(-1, 1)\n",
    "            X = np.concatenate((ones, X), axis = 1)\n",
    "        self.X = X\n",
    "        self.N, self.D = X.shape\n",
    "        self.y = y\n",
    "        self.K = len(np.unique(y))\n",
    "        self.n_iter = n_iter\n",
    "        self.lr = lr\n",
    "        \n",
    "        ### Fit B ###\n",
    "        B = np.random.randn(self.D*self.K).reshape((self.D, self.K))\n",
    "        self.I = make_I_matrix(self.y)\n",
    "        for i in range(n_iter):\n",
    "            Z = np.dot(self.X, B)\n",
    "            P = softmax_byrow(Z)\n",
    "            gradient = np.dot(self.X.T, self.I - P)\n",
    "            B += lr*gradient\n",
    "        \n",
    "        ### Return Values ###\n",
    "        self.B = B\n",
    "        self.Z = np.dot(self.X, B)\n",
    "        self.P = softmax_byrow(self.Z)\n",
    "        self.yhat = self.P.argmax(1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multiclass model is instantiated and fit below. The `yhat` value returns the class with the greatest estimated probability. We are again able to classify all observations correctly.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In-sample accuracy: 1.0\n"
     ]
    }
   ],
   "source": [
    "multiclass_model = MulticlassLogisticRegression()\n",
    "multiclass_model.fit(X, y, 10**4, 0.0001)\n",
    "print('In-sample accuracy: '  + str(np.mean(multiclass_model.yhat == y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show the distribution of our estimates of the probability that each observation belongs to the class it actually belongs to. E.g. for observations of class 1, we plot $P(y_n = 1)$. The fact that most counts are close to 1 shows that again our model is confident in its predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+MAAAFXCAYAAADJQJ+aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debwkdX3v/9dbBtwNiwckjASMoKCJqBMuCVERRHGJYNxjzJCQEKPxyjWJosnNT01uokkUydXfjQjKxA0QNXBdUBxBoyI6CAg4KotER0Zm2FRMlMXP/aPqSE/TZ87WXd3nzOv5ePTjdFV9q+tT1dWfPp/qb1WlqpAkSZIkSd25x7gDkCRJkiRpW2MxLkmSJElSxyzGJUmSJEnqmMW4JEmSJEkdsxiXJEmSJKljFuOSJEmSJHVsUcV4kqOTVJKHDpi2op32ugHt95rnMv5gMXEuZUkelOTsJDe12+64ES9vuyR/kuTLSX6U5NYkX0ny0iTb9bU9pI3pSaOMaRySHJfktweMf12SibsfYJJTk1w75NdMkhclWZvkxiS3J9mQ5LQkTxzlshcrySOSfKrdf29M8u4kO487rmnmztEzd46HuXPp5s4kj0zyjiQXJbltQt8vc+eImTvHw9y5pHPnHyX5eJLvJflxksuT/EWSHeYyf9e/jH8M+HVg4zzmORrYZpMi8NfAE4BjaLbdaaNaUJLtgbOBE4HPA88GngV8DjgBOCvJilEtf8IcB9wtKQIn07wPk+ZvaN6roWi/AM8A1gDX0ux/hwGvBu4FrE3yC8Na3jAl+UXgfODewHOAlwFPAj6aZKn2BjJ3zp+5czzMnUs0dwKPBZ4GfAdYN+ZYhsXcOX/mzvEwdy7d3PnXwPeBVwDPAE6n2T7vm8vMne7gVbUZ2NzlMhcjyT2r6qdjDmM/4NKq+sgwXmyWdfpLmi/io6rqrJ7x5yb5HPBvbZvXDyOWYUgSYPuquq2L5VXVBmBDF8uaj6q6esgv+RqaQvY5VfWhvmnvS/Jk4PYhL3NY/gLYHvitqroFIMl1wGeBo4APjzG2BTF3Loi5cyvMnQ1z5xbeU1VrAJL8LZNZAMyLuXNBzJ1bYe5smDu38Jg210w7r91PXp/kIVV1zVbnrqoFP2iOHhbw0AHTVrTTXjeg/V49434HuBi4FfgBcBnwx+2089v2vY/ze+Y9EPh0O++PgbXAgQNieQXNUZafAF8GfqMdPnVAbI8HPgjcAlzSTvs14EyaD8N/Ad8E/g64d99yzqc5sncEcEnb9mLgv7Xb4+9ojs7eBJwK3Hcr23avAev+8203l3Vvl7GB5gv1i208J86wvHu26/yxrcT0ceBm4J7t8CFtTM9ul3Uz8EOaI0G7DHgP1rcx3Exz1P1ZfW1+G/gS8J9tLB8E9uxrcy3wXpqj1t+g+WA+v92mbx4Q8/PbGA+Y63vZLqN/u5/aTnsdUH3LeADwNuA64Kfta/4PID1tprfVM9u2N9D8g/BeYMf5bqsB63kqcO2A/eePgTfQ7He3AP8XWDnLa+3QLvejc8wDWyy7Hfd64Ks0n+kbgM8AB/W1uR/wv2l+hfkpcD3NPv3wRW6Lq4H3DRj/H8CaxeS8YT0wd5o7zZ1g7txi2e24seXOvmX8bf/7NQkPzJ3mTnMnmDu3WHY7biJyZ8/rPLXdHo+bre2wfhnfbkA3ku0GtuyR5Ddpdop/pvk16x7Aw4Ed2yYvbadvR/PmQvOhI8mv0vzS9XXuSmjHA59NclBVXdq2+0PgrcApNB+yXwbe37OMfu8DPkBzdGZ6nfakSXKnAj8CHkHTJeEhwAv65n8o8I/A/6JJWP9A0wXn7Pb1jqY56viPwCbgVTPEsZEmmb0DuLPdFgAb57rurV+g6WL0T8BraXauQR7btj17hum0054KPAa4oGf8W2l25hcC+9AkmV8EngiQ5EXAm2k+nP9O0334V4Gfn8Ob5CXA/wHe3ba7P00C+mySX62qH/Us74nAATQfvE00SewM4HeSvKqq7uxp+7vA5VV1STs8l/fyWTRfAJe2McAMR9bbbs8fa7fJX9N8qT8deAswRbPNe50IfJTmn4GH0ewfdwKr57qt5uk1NF+IfwDs2r72+2i6oM1kFc3nY2v7wmz2oOlitgG4L8378Lkkq6rqa22bE2i+JF4LXAnsAhzcLntB2yLJvYG9abp19bsC2H8R6zQK5s67mDvNnebOMeXOJcjceRdzp7nT3Dl5ufMJwM+Ab83acr6Vfl/VfzSDj6L1PmY8Qgn8OXDTLMs4H/j8gPFn0hxx2bFn3ANojlR9uB2+B/Bd4OMDjoT9/KhTX2wnzBJPaJLb77YbeZe+WG8HHtIz7pnt636673U+DHx7Dtv48/QclZ3ruvccOSrgyDksZ/pI3lO20uaIts3z2uFD2uFz+tq9qB1/WDv8NuCrW3nd+9EcyXpX3/i9gNuA43rGXUtzBPNBfW0P7o+fJindDrxqAe/ltcB7B8zzOnqOUNKcG1LA0X3tTqY56vbAvm21pq/d22iOnGcu22or2/BUBh+h/Gxfuz9vx//iYvaFrS17wPTt2u38TXqOkAOXA2/Zynzz3hY0X8YFvGTAtPcCV893247igbnT3GnuNHdOUO4c8BqT/su4ufOuWM2d5k4wd/ZOH1vubF/nV2kOQr1zLu2HdTGjZ9F0w+h9HDSH+b4C7JTkvUmekWSmo4aDPJ6mO8Mt0yOq6oc0R1We0I5a2T4+2DfvWcAdM7zu3c6RSfKAJG9KcjXNjn478B6aD9U+fc2/VVueG/CN9u8n+9p9A1jZnlMwX3NZ92l30BwRm81c4pipzRl9wx+kSTLT55t9BTggyf9O8qQk9+lr/+s0Sf19aa6GuqI94r2BZjs9vq/9l6rq+70jquoLNN2TX9wz+gU0X4w/v4DCPN/LuXh8u64f6Bv/XppuN/3n3H2sb/gymq5au7XDs22r+Rq0PGiO1I5MG/t5SW6k2QdvB/alOSo77SvA0Ulem2RV+q6aysK2xfQ+WluZNknMnXcxd5o7wdw5rty51Jg772LuNHeCuXMicmeS3Wk+71cDr5zLPMMqxi+vqnW9D+Ci2Waqqs8CzwUeTJOMNif5dNsdZjY7M/jqmN8Hdmqf797+3dS33DtpzicYZNBrvht4CU23psNpkv7L2mn36mt7c9/wbVsZv4I5dKsaYC7rPm1Tbdl9Zibfbf/utZU2v9TXdtr1vQPVXNTiZpouIwD/CvwJzTlMnwRuSvLh3HWrkV3bv5+m+fD0Pn6FphtJr5muivpe4FlJ7tcOvxj4TFV9r6fNfN7LudiZ5ih7/8VJvt8zvddNfcPT800ve7ZtNV+zLW+Q6ff3l7bSZkZJHkPT3epWmqthHkSznS/tW+7LabrD/QFNAtyU5ISe5LeQbXEzTSE+qEvRTtx9e4ybufMu5k5zJ5g7x5U7lxpz513MneZOMHeOPXcm2QU4l+ZAy1Nqy1MdZjT22/xU1ZlV9QSaD/OzaBLZOZn9FkQ3AQ8aMP5B3LUjTH94du1t0B4JeeBMIfW1vRdwJPCPVXViVX22TfoznQPThbms+7RBvxAOso7mvKhnbqXNM2m69Xy1b/xuvQNp7qu3E/A9aPrWVNU7qupAmu2+muZCIKe3s9zY/j2aux/p/jXg2Dmu03uA+9Akxn3bed/TE9co3subgJ1z93sJTr8/NzIPc9hWXVhH0x3ttxY4/7Npjkr+dlX9W1Vd2G7nLb6wq+rWqnpNVT2U5sv474A/Bf6/dvq8t0VV/SdNV69HDJi8P835bsuCuXNBzJ2DmTuHY8nmzm2JuXNBzJ2DmTuHY8nnziQPoCngdwGe1HdAZqvGXoxPazfQR2mOWOzOXUelfkpzAn2/zwJPT3L/6RHt899qp0HT3WQDzVHQXkcx99u63ZPmKGL/5fSPnuP8ozCXdZ+X9gjbPwNPS3Jk//R23FNpzr3oPxr3vL7h59LsWxf0jaeqbq6q02m6GD2yHf1FmotaPLT/SHf7+OYc1+Hqdpkvbh8/ZsvbWM3nvZxpv+v3WZp17d/HXkRzFPpLc3iNgWbYViPXHmF+M/CMJM8e1CbJ4VvpvnMfmouDVE/7Q9lKF6Wq+o+qejNNd6a7res8t8XZNJ+Pn9+PMs1Fe36JxV0cZCKZO+fF3Dl4HcydQ7AMcuc2xdw5L+bOwetg7hyCpZ4727g+RnMB4SdX1VVba9+v0/uM90vyBpqjW+fRXJ5/JfDfaW7tMH0Vwa8DL03yfJr+9z9qPyR/Q3MRg7VJ3kTzBrya5g15A0BV/SzJ64F3JjmZ5pySh9BcAfIHNOdcbFVV/SDJl4A/S7KRppvRH3BXV5hxmHXdF+gNNFc0PCPJ24FPtK99BE3XjnNoLurS7xFJ3k1z9cx9aa7o+dmqWguQ5CSapHcBTdetfWmS1qegOe8oyV8Ab08y1S73BzTb+Ak0FxJ5/xzX4V+Bt9N0M/pIVd06PWGe7+XXgccleQZN158bquraAe0+QXOxk39pY7+C5p6Zfwj8fVXN1C1toNm2VYf+HngUcHqSU2luTXETzWf02TQXo+nvmjbtHOA44NR2v9gX+J+0R6ynJbmApji+jKZr0RPaZU7f53ah2+IfaS6OcnaSv6e5Wus/0NxeZij3TR03c+eCmTtnZu4cjiWbO9t/KJ/WDj68Hfecdvja9pemJc3cuWDmzpmZO4djyeZO4EM0F/R7BXDfJL3Xr7i6trwH+d3V4q4WdzQs/H6PNJfi/yRNt56f0pwzcAo9V9yj6Xbx8XbjFFve7/G/Mbf7PR5Hc4/hn9B0hfhNmnNLTpjjuuxFs/P/qH1z3tbGXsAhPe3Op+8KnNx1dcE/7Bv/unb8ilm28d2uajnXdae93+M839MVNOeyfKV93R+32+xP+2Plris1/na7rFvabfR+2qs5tu1Wt9tmU/s+f5vm9gIP6Hu9p9F8Qf6QpgvPVcC7gP172lzLgKtN9kzfqV1G0RydWuh7+XCa2xr8Jz1XQKXvqpbtuOn7PW6kOSr5LWa+3+OTZvgM7TWfbTVgvU5l8FUt+/e76TgO2drrtW1DU9SeR/N5uZ3miP8H6LlvYv+y23Evb2P/r3ZfelK7Xr2f3zfR3A/1B+1+dhnw3+e738wQ+6/QnLfz4zb2U+m7B+k4H5g7zZ3mTjB3brHsdtzYcmfPug96nDqfz8OoHpg7zZ3mTjB3brHsdtw4c+dMebPou+r9oMf0Ze23KUl+jeZXst+rqvfM1l6SZO6UpIUwd0qaybIvxpPsTXPE7d9pjnztR3Oz99uAR1ZzwSdJUg9zpyTNn7lT0nyM9ZzxjvwXzYn3v0fTleRmmm42x5sQJWlG5k5Jmj9zp6Q5W/a/jEuSJEmSNGkm5tZmkiRJkiRtKyzGJUmSJEnqWKfnjB9xxBF1zjnndLlISduGjDuAUTJ3ShqRZZs7zZuSRmSoebPTX8ZvuGFe96GXJGHulKT5Mm9KWgrspi5JkiRJUscsxiVJkiRJ6pjFuCRJkiRJHbMYlyRJkiSpYxbjkiRJkiR1zGJckiRJkqSOWYxLkiRJktQxi3FJGpEkOyY5M8k3kqxP8utJdk5ybpIr2787jTtOSZIkdc9iXJJG50TgnKp6OPAoYD1wPLC2qvYB1rbDkiRJ2sZYjEvSCCR5APB44BSAqrqtqm4BjgTWtM3WAEeNJ0JJkiSNk8W4JI3GQ4DNwLuTXJzk5CT3BXarqo0A7d9dxxmkJEmSxmPFuAOQtHy8891XcMONty3qNR64yw780e8/YkgRjdUK4DHAy6vqwiQnMo8u6UmOBY4F2HPPPUcToaSJYO6UpPlbDrlzTsV4kmuBHwF3AndU1aokOwOnA3sB1wLPq6qbRxOmpKXghhtv4+Cn7reo1/jCJ9YPKZqx2wBsqKoL2+EzaYrx65PsXlUbk+wObBo0c1WdBJwEsGrVquoiYEnjYe6UpPlbDrlzPt3Un1hVB1TVqnbYixBJ0gyq6vvAd5M8rB11GPB14GxgdTtuNXDWGMKTJEnSmC2mm/qRwCHt8zXA+cCrFxmPJC0nLwfel2QH4Brg92kOgp6R5BjgO8BzxxifJEmSxmSuxXgBn0pSwDva7pNbXIQoiRchkqQeVXUJsGrApMO6jkWSJEmTZa7F+MFVdV1bcJ+b5BtzXYAXIZIkSZIkaUtzOme8qq5r/24CPgIcSHsRIoDZLkJUVauqatXU1NRwopYkSZIkaQmbtRhPct8k959+DjwZuBwvQiRJkiRJ0oLMpZv6bsBHkky3f39VnZPkK3gRIkmSJEmS5m3WYryqrgEeNWD8jXgRIkmSJEmS5m0+9xmXJEmSJElDYDEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSdLESLJdkouTfLQd3jvJhUmuTHJ6kh3GHaMkDYPFuCRJkibJK4D1PcNvAk6oqn2Am4FjxhKVJA2ZxbgkSZImQpKVwNOBk9vhAIcCZ7ZN1gBHjSc6SRoui3FJkiRNircCrwJ+1g7vAtxSVXe0wxuAPcYRmCQNm8W4JEmSxi7JM4BNVXVR7+gBTWuG+Y9Nsi7Jus2bN48kRkkaJotxSZIkTYKDgWcmuRY4jaZ7+luBHZOsaNusBK4bNHNVnVRVq6pq1dTUVBfxStKiWIxLkiRp7KrqNVW1sqr2Al4AfKaqXgScBzynbbYaOGtMIUrSUFmMS5IkaZK9GnhlkqtoziE/ZczxSNJQrJi9iSRJktSdqjofOL99fg1w4DjjkaRR8JdxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdcxiXJIkSZKkjlmMS5IkSZLUMYtxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdcxiXJIkSZKkjlmMS5IkSZLUsRXjDkCSlqsk1wI/Au4E7qiqVUl2Bk4H9gKuBZ5XVTePK0ZJkiSNh7+MS9JoPbGqDqiqVe3w8cDaqtoHWNsOS5IkaRtjMS5J3ToSWNM+XwMcNcZYJEmSNCYW45I0OgV8KslFSY5tx+1WVRsB2r+7ji06SZIkjY3njEvS6BxcVdcl2RU4N8k35jpjW7wfC7DnnnuOKj5JkiSNib+MS9KIVNV17d9NwEeAA4Hrk+wO0P7dNMO8J1XVqqpaNTU11VXIkiRJ6ojFuCSNQJL7Jrn/9HPgycDlwNnA6rbZauCs8UQoSZKkcbKbuiSNxm7AR5JAk2vfX1XnJPkKcEaSY4DvAM8dY4ySJEkaE4txSRqBqroGeNSA8TcCh3UfkSRJkiaJ3dQlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdWzOxXiS7ZJcnOSj7fDeSS5McmWS05PsMLowJUmSJElaPubzy/grgPU9w28CTqiqfYCbgWOGGZgkSZIkScvVnIrxJCuBpwMnt8MBDgXObJusAY4aRYCSJEmSJC03c/1l/K3Aq4CftcO7ALdU1R3t8AZgjyHHJkmSJEnSsjRrMZ7kGcCmqrqod/SApjXD/McmWZdk3ebNmxcYpiRJkiRJy8dcfhk/GHhmkmuB02i6p78V2DHJirbNSuC6QTNX1UlVtaqqVk1NTQ0hZEmSJEmSlrZZi/Gqek1VrayqvYAXAJ+pqhcB5wHPaZutBs4aWZSSJEmSJC0ji7nP+KuBVya5iuYc8lOGE5IkSZIkScvbitmb3KWqzgfOb59fAxw4/JAkSZIkSVreFvPLuCRJkiRJWgCLcUmSJEmSOmYxLkmSJElSxyzGJUmSJEnqmMW4JEmSJEkdsxiXJEmSJKljFuOSJEmSJHXMYlySJEmSpI5ZjEuSJEmS1DGLcUmSJEmSOmYxLkmSJElSxyzGJUmSJEnqmMW4JEmSJEkdsxiXJEnSREhyryRfTnJpkiuSvL4dv3eSC5NcmeT0JDuMO1ZJWiyLcUmSJE2KnwKHVtWjgAOAI5IcBLwJOKGq9gFuBo4ZY4ySNBQW45IkSZoI1bi1Hdy+fRRwKHBmO34NcNQYwpOkobIYlyRJ0sRIsl2SS4BNwLnA1cAtVXVH22QDsMe44pOkYbEYlyRJ0sSoqjur6gBgJXAgsN+gZv0jkhybZF2SdZs3bx51mJK0aBbjkiRJmjhVdQtwPnAQsGOSFe2klcB1A9qfVFWrqmrV1NRUd4FK0gJZjEuSJGkiJJlKsmP7/N7Ak4D1wHnAc9pmq4GzxhOhJA3PitmbSJIkSZ3YHViTZDuaH43OqKqPJvk6cFqSvwUuBk4ZZ5CSNAwW45I0Qu0/lOuA71XVM5LsDZwG7Ax8FXhxVd02zhglaVJU1deARw8Yfw3N+eOStGzYTV2SRusVNF0sp3mvXEmSJFmMS9KoJFkJPB04uR0O3itXkiRJWIxL0ii9FXgV8LN2eBe8V64kSZKwGJekkUjyDGBTVV3UO3pA07vdK7ed3/vlSpIkLWMW45I0GgcDz0xyLc0F2w6l+aV81nvlgvfLlSRJWu4sxiVpBKrqNVW1sqr2Al4AfKaqXoT3ypUkSRIW45LUtVcDr0xyFc055N4rV5IkaRvkfcYlacSq6nzg/Pa598qVJEmSv4xLkiRJktQ1i3FJkiRJkjpmMS5JkiRJUscsxiVJkiRJ6pjFuCRJkiRJHbMYlyRJkiSpYxbjkiRJkiR1zGJckiRJkqSOWYxLkiRJktQxi3FJkiRJkjpmMS5JkiRJUscsxiVJkiRJ6tisxXiSeyX5cpJLk1yR5PXt+L2TXJjkyiSnJ9lh9OFKkiRJkrT0zeWX8Z8Ch1bVo4ADgCOSHAS8CTihqvYBbgaOGV2YkiRJkiQtH7MW49W4tR3cvn0UcChwZjt+DXDUSCKUJEmSJGmZmdM540m2S3IJsAk4F7gauKWq7mibbAD2GE2IkiRJkiQtL3Mqxqvqzqo6AFgJHAjsN6jZoHmTHJtkXZJ1mzdvXnikkiRJkiQtE/O6mnpV3QKcDxwE7JhkRTtpJXDdDPOcVFWrqmrV1NTUYmKVJEmSJGlZmMvV1KeS7Ng+vzfwJGA9cB7wnLbZauCsUQUpSZIkSdJysmL2JuwOrEmyHU3xfkZVfTTJ14HTkvwtcDFwygjjlCRJkiRp2Zi1GK+qrwGPHjD+GprzxyVJkiRJ0jzM65xxSZIkSZK0eBbjkiRJkiR1zGJckiRJkqSOWYxLkiRJktQxi3FJkiRJkjpmMS5JkiRJUscsxiVJkiRJ6pjFuCRJkiRJHbMYlyRJkiSpYxbjkiRJkiR1zGJckiRJkqSOWYxLkiRJktQxi3FJkiRJkjpmMS5JkiRJUscsxiVpBJLcK8mXk1ya5Iokr2/H753kwiRXJjk9yQ7jjlWSJEndsxiXpNH4KXBoVT0KOAA4IslBwJuAE6pqH+Bm4JgxxihJkqQxsRiXpBGoxq3t4Pbto4BDgTPb8WuAo8YQniRJksbMYlySRiTJdkkuATYB5wJXA7dU1R1tkw3AHjPMe2ySdUnWbd68uZuAJUmS1BmLcUkakaq6s6oOAFYCBwL7DWo2w7wnVdWqqlo1NTU1yjAlSZI0BhbjkjRiVXULcD5wELBjkhXtpJXAdeOKS5IkSeNjMS5JI5BkKsmO7fN7A08C1gPnAc9pm60GzhpPhJIkSRqnFbM3kSQtwO7AmiTb0Rz4PKOqPprk68BpSf4WuBg4ZZxBSpIkaTwsxiVpBKrqa8CjB4y/hub8cUmSJG3D7KYuSZKkiZDkwUnOS7I+yRVJXtGO3znJuUmubP/uNO5YJWmxLMYlSZI0Ke4A/qyq9qO56OXLkuwPHA+srap9gLXtsCQtaRbjkiRJmghVtbGqvto+/xHNhS/3AI4E1rTN1gBHjSdCSRoei3FJkiRNnCR70Vx740Jgt6raCE3BDuw6oP2xSdYlWbd58+YuQ5WkBbEYlyRJ0kRJcj/gQ8BxVfXDucxTVSdV1aqqWjU1NTXaACVpCCzGJUmSNDGSbE9TiL+vqj7cjr4+ye7t9N2BTeOKT5KGxWJckiRJEyFJgFOA9VX1lp5JZwOr2+ergbO6jk2Shs37jEuSJGlSHAy8GLgsySXtuNcCbwTOSHIM8B3guWOKT5KGxmJckiRJE6GqPg9khsmHdRmLJI2a3dQlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdcxiXJIkSZKkjlmMS5IkSZLUMYtxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSerYrMV4kgcnOS/J+iRXJHlFO37nJOcmubL9u9Pow5UkSZIkaembyy/jdwB/VlX7AQcBL0uyP3A8sLaq9gHWtsOSJEmSJGkWsxbjVbWxqr7aPv8RsB7YAzgSWNM2WwMcNaogJUmSJElaTuZ1zniSvYBHAxcCu1XVRmgKdmDXYQcnSZIkSdJyNOdiPMn9gA8Bx1XVD+cx37FJ1iVZt3nz5oXEKEmSJEnSsjKnYjzJ9jSF+Puq6sPt6OuT7N5O3x3YNGjeqjqpqlZV1aqpqalhxCxJkiRJ0pI2l6upBzgFWF9Vb+mZdDawun2+Gjhr+OFJkiRJkrT8rJhDm4OBFwOXJbmkHfda4I3AGUmOAb4DPHc0IUqSJEmStLzMWoxX1eeBzDD5sOGGI0mSJEnS8jevq6lLkiRJkqTFsxiXJEmSJKljFuOSJEmSJHXMYlySJEmSpI5ZjEuSJEmS1DGLcUmSJEmSOmYxLkmSJElSxyzGJWkEkjw4yXlJ1ie5Iskr2vE7Jzk3yZXt353GHaskSZK6ZzEuSaNxB/BnVbUfcBDwsiT7A8cDa6tqH2BtOyxJkqRtjMW4JI1AVW2sqq+2z38ErAf2AI4E1rTN1gBHjSdCSZIkjZPFuCSNWJK9gEcDFwK7VdVGaAp2YNfxRSZJkqRxsRiXpBFKcj/gQ8BxVfXDecx3bJJ1SdZt3rx5dAFKkiRpLCzGJWlEkmxPU4i/r6o+3I6+Psnu7fTdgU2D5q2qk6pqVVWtmpqa6iZgSZIkdcZiXJJGIEmAU4D1VfWWnklnA6vb56uBs7qOTZIkSeO3YtwBSNIydTDwYuCyJJe0414LvBE4I8kxwHeA544pPkmSJI2RxbgkjUBVfR7IDJMP6zIWSZIkTR67qZNtgpEAABNISURBVEuSJEmS1DGLcUmSJEmSOmYxLkmSJElSxyzGJUmSJEnqmMW4JEmSJEkdsxiXJEmSJKljFuOSJEmSJHXMYlySJEmSpI5ZjEuSJEmS1DGLcUmSJEmSOmYxLkmSJElSxyzGJUmSJEnqmMW4JEmSJEkdsxiXJEmSJKljFuOSJEmSJHXMYlySJEmSpI5ZjEuSJEmS1DGLcUmSJE2EJO9KsinJ5T3jdk5ybpIr2787jTNGSRoWi3FJkiRNilOBI/rGHQ+srap9gLXtsCQteRbjkiRJmghV9Tngpr7RRwJr2udrgKM6DUqSRsRiXJIkSZNst6raCND+3XXM8UjSUFiMS5IkaclLcmySdUnWbd68edzhSNKsLMYlSZI0ya5PsjtA+3fToEZVdVJVraqqVVNTU50GKEkLYTEuSZKkSXY2sLp9vho4a4yxSNLQWIxLkiRpIiT5AHAB8LAkG5IcA7wRODzJlcDh7bAkLXkrxh2AJEmSBFBVL5xh0mGdBiJJHZj1l/Ek70qyKcnlPeN2TnJukivbvzuNNkxJkiRJkpaPuXRTPxU4om/c8cDaqtoHWNsOS5IkSZKkOZi1GK+qzwE39Y0+EljTPl8DHDXkuCRJkiRJWrYWegG33apqI0D7d9fhhSRJkiRJ0vI28qupJzk2ybok6zZv3jzqxUmSJEmSNPEWWoxfn2R3gPbvppkaVtVJVbWqqlZNTU0tcHGSJEmSJC0fCy3GzwZWt89XA2cNJxxJkiRJkpa/udza7APABcDDkmxIcgzwRuDwJFcCh7fDkiRJkiRpDlbM1qCqXjjDpMOGHIskLRtJ3gU8A9hUVY9sx+0MnA7sBVwLPK+qbh5XjJIkSRqfkV/ATZK2UacCR/SNOx5YW1X7AGvbYUmSJG2DLMYlaQSq6nPATX2jjwTWtM/XAEd1GpQkSZImhsW4JHVnt6raCND+3XXM8UiSJGlMLMYlaQIlOTbJuiTrNm/ePO5wJEmSNGQW45LUneuT7A7Q/t00U8OqOqmqVlXVqqmpqc4ClCRJUjcsxiWpO2cDq9vnq4GzxhiLJEmSxshiXJJGIMkHgAuAhyXZkOQY4I3A4UmuBA5vhyVJkrQNmvU+45Kk+auqF84w6bBOA5EkSdJE8pdxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdcxiXJIkSZKkjlmMS5IkSZLUMYtxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdcxiXJIkSZKkjlmMS5IkSZLUMYtxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMrxh2AJEnSOLzz3Vdww423Leo1HrjLDvzR7z9iSBFJkrYlFuOSJGmbdMONt3HwU/db1Gt84RPrhxSNJGlbYzd1SZIkSZI6ZjEuSZIkSVLH7KYuSZK0QF+9eBN//08XL+o1Lv7a5kV3l5ckLT0W45IkSQv0k9t+tuhC+osXfn9I0UiSlhK7qUuSJEmS1LGJ/2V8sbcd8ZYj2hZ4ex5JkiRpaZn4Ynyxtx3xliPaFnh7Hklz5UFuSZImw8QX45IkaXg8yC1J0mTwnHFJkiRJkjrmL+OSJEna5g3j+ivf/NaNPGzfXRb1Gp4KMlxeV2dL7ueTZdkX48O4/6c7y2RaTslkseviPWrVy388tC0Yxn5u7lSvYVx/5YsXnuc1XCaM19XZkvv5ZFlUMZ7kCOBEYDvg5Kp641CiGqJh3P/TnWUyLadksth18R61S8uoc6f/eGhbMJzvAHPnUrIU/u+UpPlY8DnjSbYD3g48FdgfeGGS/YcVmCQtR+ZOSZo/c6ek5WgxF3A7ELiqqq6pqtuA04AjhxOWJC1b5k5Jmj9zp6RlZzHF+B7Ad3uGN7TjJEkzM3dK0vyZOyUtO6mqhc2YPBd4SlX9YTv8YuDAqnp5X7tjgWPbwYcB31x4uIv2QOCGMS5/MYx9PJZq7Es1blhY7DdU1RGjCGbYJjB3Tsq+MilxgLHMZFJimZQ4YOnHsqxy54T9z7kYk7RfDZvrtjS5bncZat5czAXcNgAP7hleCVzX36iqTgJOWsRyhibJuqpaNe44FsLYx2Opxr5U44alHfscTVTunJTtPSlxgLHMZFJimZQ4wFg6NmvunKT/ORdjOb+XrtvS5LqNzmK6qX8F2CfJ3kl2AF4AnD2csCRp2TJ3StL8mTslLTsL/mW8qu5I8qfAJ2luMfGuqrpiaJFJ0jJk7pSk+TN3SlqOFnWf8ar6OPDxIcXShaXcdcnYx2Opxr5U44alHfucTFjunJTtPSlxgLHMZFJimZQ4wFg6NWG5c5SW83vpui1NrtuILPgCbpIkSZIkaWEWc864JEmSJElagGVTjCc5Isk3k1yV5PgB038pydokX0tyfpKVPdNWJ7myfazuNvKFx57kgCQXJLminfb8pRB3z/QHJPlekrd1F/XPl72Y/WXPJJ9Ksj7J15PstYRi/4d2f1mf5J+TpMO435VkU5LLZ5ieNqar2tgf0zNtrJ/RpWi2/aRt87x2H74iyft7xt+Z5JL2segLJM1hnz2hZ3nfSnJLz7ShvveLjKXr7bJnkvOSXNx+Jp7WM+017XzfTPKUccSRZK8k/9WzTf5lMXHMMZbOvssXGcvQ9hVz5/IxKZ/5UZikPDJsk5SXhm1S8twoLJncWVVL/kFzIY+rgYcAOwCXAvv3tfkgsLp9fijwnvb5zsA17d+d2uc7LZHY9wX2aZ//IrAR2HHS4+6ZfiLwfuBtS2V/aYfPBw5vn98PuM9SiB34DeAL7WtsB1wAHNJh7I8HHgNcPsP0pwGfAAIcBFzYjh/rZ3QpPua4n+wDXDy9LYFde6bd2mUsfe1fTnNhpqG/94uJZRzbheY8tj9pn+8PXNvz/FLgnsDe7etsN4Y49prp8zzCWDr5Ll9MLCPYV8ydy+AxKZ/5CVy3oeaRMa3bRNYYo1y3dnhoeW5E67ckcudy+WX8QOCqqrqmqm4DTgOO7GuzP7C2fX5ez/SnAOdW1U1VdTNwLjC0G7nPwYJjr6pvVdWV7fPrgE3AVCdRL26bk+SxwG7ApzqItd+CY0+yP7Ciqs4FqKpbq+o/uwkbWNx2L+BeNAn3nsD2wPUjj3h64VWfA27aSpMjgX+txpeAHZPszvg/o0vRXPaTPwLe3m5TqmrTGGPp9ULgA+3zYb/3i4ll2OYSSwEPaJ//AnfdU/lI4LSq+mlVfRu4qn29ruMYtkn6Ll/Ud9wwmTuXjUn5zI/CJOWRYZukvDRsE5PnRmGp5M7lUozvAXy3Z3hDO67XpcCz2+fPAu6fZJc5zjtKi4n955IcSFNkXT2iOPstOO4k9wDeDPzFyKMcbDHbfF/gliQfbrta/WOS7UYe8V0WHHtVXUCTSDe2j09W1foRxzsfM63buD+jS9Fcttm+wL5JvpDkS0l6v2julWRdO/6oDmIBmu5wNL/6fGa+83YQC3S/XV4H/G6SDTRXkH75PObtIg6Avdtc+Nkkj1tgDPOJpavv8sV+Nw9zX5mNuXNpmJTP/ChMUh4ZtknKS8O2lPLcKExE7lwuxfig8177LxP/58ATklwMPAH4HnDHHOcdpcXE3rxAcxTnPcDvV9XPRhVon8XE/VLg41X1XcZjMbGvAB7XTv81mq49R48s0rtbcOxJHgrsB6ykSSqHJnn8KIOdp5nWbdyf0aVoLttsBU1X9UNofgE+OcmO7bQ9q2oV8DvAW5P88ohjmfYC4MyqunMB8446Fuh+u7wQOLWqVtJ0p3tPezBzmNtlMXFspNkmjwZeCbw/yQNYuEn6Ll/sd/Mw95XZmDuXhkn5zI/CJOWRYZukvDRsSynPjcJE5M7lUoxvAB7cM7ySvu4vVXVdVf12+2H/y3bcD+Yy74gtJnbahPUx4K/aLhZdWUzcvw78aZJrgX8Cfi/JGzuJurHY/eXitkvPHcC/0ZyP0pXFxP4s4Ett1/pbac6TOaibsOdkpnUb92d0KZrLNtsAnFVVt7fdHr9JU5xPn/ZCVV1Dc42ER484lmkvYMtu4cN+7xcTyzi2yzHAGe0yL6A5zeSBc5x35HG0XWZvbMdfRNMza98FxjGnWDr8Ll/Ud/OQ95WFxmrunCyT8pkfhUnKI8M2SXlp2JZSnhuFycidNQEn2C/2QfMLzzU0XQqnL0DwiL42DwTu0T7/X8Ab6q6T9L9Nc4L+Tu3znZdI7DvQnMdx3FLa5n1tjqb7C7gtZptv17afaoffDbxsicT+fODT7Wts3+47v9Xxtt+LmS+k8XS2vJDGl9vxY/2MLsXHHPeTI4A1PfvMd4Fd2m18z57xV7KVi5wNI5a23cOAa4H0jBvqe7/IWDrfLu3n4ej2+X40/wwEeARbXszpGhZ+AbfFxDE1vVyaXkLfG/X7Q0ff5YuMZaj7Svs6e2HuXNKPSfnMT+C6DTWPjGndJrLGGPG6DT3PjWgdJz53jn0jDXFjPw34Fs0Rtb9sx70BeGb7/DntjvIt4OTpHaid9gc0F8O4iqar95KIHfhd4Hbgkp7HAZMed99rHE3HxfgQ9pfDga8BlwGnAjsshdhpDiS8A1gPfB14S8dxf4CmO9rtNEcdjwFeAryknR7g7e16XQas6pl3rJ/RpfiYw34S4C3tvnAZ8IJ2/G+0w5e2f48ZdSzt8OuANw6Yd6jv/UJjGcd2oblwzhfaZV4CPLln3r9s5/sm8NRxxEFzHuEV7fivMoSDe3OIpbPv8oXGMux9BXPnsnlMymd+ktZtFHlkDOs2sTXGqNZt2HluROu2JHJn2gVKkiRJkqSOLJdzxiVJkiRJWjIsxiVJkiRJ6pjFuCRJkiRJHbMYlyRJkiSpYxbjkiRJkiR1zGJckiRJkqSOWYxLkiRJktQxi3GNRJKdklyf5JfHHUu/JGcmeeW445CkfuZOSZo/c6eWKotxLUiSzySp9nFHkquT/HFPk9cCH6+qq8cQ20uTfDvJT5JclORxfU1eD/xVkl/oOjZJ27ZJzZ1JHp/k7CTfa2M7ekAzc6eksZjg3PmaJF9J8sMkm5P83ySP7Gtm7tSMLMa1UI8GXgfsDuwDfAL4P0keneQ+wB8Cp3QdVJLnAycCf9fG+EXgE0n2nG5TVZcB1wC/23V8krZ5E5k7gfsBlwOvAP5rUANzp6QxmtTceQjw/wO/ARwK3AF8OsnO0w3Mndoai3HNW9sFaEfgC1X1/ar6NvA3QGiS5dOAnwFf6JlnQ38XnSS/0v56vf8Qw3slcGpVvbOq1lfVy4GNwJ/0tTsbeOEQlytJWzXJubOqPl5Vr62qM9sYZmLulNSpCc+dT6mqd1fV5W3R/WJgCji4r6m5UwNZjGshHtv+vbRn3Mr27ybgccBFVVU90y8Afq3vdd4KnFxVX+8dmeS1SW6d5dHf9ZwkO7Sxfapv0qdojlj2+jJwYJJ7z7KukjQsE5k758ncKalrSyl33p+mvrq5b7y5UwOtGHcAWpIeC1xXVZvh50csTwSuBj5N01VoY988FwAvnR5IchTN0cznDXj9fwHOmCWG7w0Y90BgO+D6vvHXA0/qG3cdsD3wi23ckjRqk5o758PcKalrSyl3nghc0i6/l7lTA1mMayEeCzwoya00xW+AjwAvrKqftEf9+gviLwFvbs+h+THwT8AbqurG/hevqpuAmxYRX/UNZ8C46XMiPUIpqSuTnjvnwtwpqWtLIncmeQvwm8BvVtWdfZPNnRrIYlwL8WjgLcA7aJLLxqrqPcfwBmCnvnkuAm4DVrXz3wG8fdCLJ3ktzVUxt+apVfXvfeNuAO4EHtQ3flfunqSnL6yxeZblSNKwTGrunA9zp6SuTXzuTHIC8ALgiVV1zYAm5k4NZDGueUmyN01C+XRVXTVDs4uBo3tHVNVPk1wM/BawGvidqrp9hvkX1F2oqm5LchFwOPDBnkmHAx/qa/5Imi5P/UW6JA3dJOfOeTJ3SurMUsidSU6kKcQPqapvzNDM3KmBLMY1X9MX0Vi3lTafBN6UZJe+7kAX0Nw259yq+uhMMy+yu9BbgPck+TLNVTVfQnN+zr/0tXsccM4ClyFJ8zXRuTPJ/YCHtoP3APZMcgBwU1V9p6epuVNSlyY9d76d5grqRwE3J5nunXlrVd3a09TcqYG8mrrm67HAtwedczOtvbXDl2mOEva6hObWE6+820xDUlWnA8cBf9Uu7zeBp1XVf0y3SXIv4FnAO0cVhyT1mejcSdOV8+L2cW/g9e3zN0w3MHdKGoNJz50vpbmC+lqai8hNP/58uoG5U1uTLe8CIA1HkiNorii5//RFLJJ8Criyql425theBhxZVU8eZxyS1M/cKUnzZ+7UUmU3dY1EVZ3Tdt1ZmeQnNOfy/Arw/LEG1rgdePm4g5CkfuZOSZo/c6eWKn8Z10glOQT4DPBN4Jiq+uJ4I5KkyWfulKT5M3dqqbEYlyRJkiSpY17ATZIkSZKkjlmMS5IkSZLUMYtxSZIkSZI6ZjEuSZIkSVLHLMYlSZIkSeqYxbgkSZIkSR2zGJckSZIkqWMW45IkSZIkdez/ATL2/IAN0wgaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1224x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize = (17, 5))\n",
    "for i, y in enumerate(np.unique(y)):\n",
    "    sns.distplot(multiclass_model.P[multiclass_model.y == y, i],\n",
    "                 hist_kws=dict(edgecolor=\"darkblue\"), \n",
    "                 color = 'cornflowerblue',\n",
    "                 bins = 15, \n",
    "                 kde = False,\n",
    "                 ax = ax[i]);\n",
    "    ax[i].set_xlabel(xlabel = fr'$P(y = {y})$', size = 14)\n",
    "    ax[i].set_title('Histogram for Observations in Class '+ str(y), size = 16)\n",
    "sns.despine()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
